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ABSTRACT: 

An  algorithm  for  rational  Chebyshev  approximation  based  on  computing 
the  zeros  of  the  error  curve  was  investigated.   At  each  iteration  the  pro- 
posed zeros  are  corrected  by  changing  them  toward  the  abscissa  of  the 
adjacent  extreme  of  largest  magnitude.   The  algorithm  is  formulated  as  a 
numerical  solution  of  a  certain  system  of  ordinary  differential  equations. 
Convergence  is  obtained  by  showing  the  system  is  asymptotically  stable  at 
the  zeros  of  the  best  approximation.   With  an  adequate  initial  guess,  the 
algorithm  has  never  failed  for  functions  which  have  a  standard  error  curve, 
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1.   Introduction 

We  will  consider  the  problem  of  computing  the  best  approximation  to  a 

m      .     n 
continuous  function   f(x)   by  a  rational  function  R   =   ^  a.x  /   y^  b .xJ 

mn    1=0  1     j=0  J 

where  the  measure  of  the  error  is  the  weighted  sup  norm, 


f  -  R    =  sup 
a<x<b 


f(x)  -  R         (x) 
mn 

w(x) 


w(x)  >  0  on  [a,b]  .   A  number  of  methods  have  been  proposed,  and  a  brief 
description  and  numerical  comparison  of  some  of  them  is  given  by  Lee  and 
Roberts  [8] . 

We  have  implemented  and  analyzed  an  algorithm  for  computing  the  best 
rational  approximation.   This  algorithm  [7]  has  not  yet  appeared  in  the 
literature.   The  convergence  rate  is  linear,  and  although  the  algorithm 
is  relatively  slow  compared  to  some  others,  the  algorithm  was  successful 
in  instances  where  others,  such  as  the  Remes  algorithm  and  Maehly's  second 
method,  fail.   For  functions  which  have  a  standard  error  curve,  the  algorithm 
has  never  failed  to  converge  to  the  best  approximation,  provided  the  initial 
approximant  did  not  have  a  pole  in  [a,b]. 

The  algorithm  is  formulated  as  the  numerical  solution  of  a  system  of 
m  +  n  +  1  differential  equations.   The  dependent  variables  are  points  of 
zero  error  in  a  proposed  approximation.   When  there  are  m  +  n  +  2   alternations 
in  the  error  curve,  the  system  of  equations  is  at  a  rest  point.   Under  the 
assumption  of  a  relatively  mild  hypothesis,  the  algorithm  is  proved  convergent 
by  showing  that  the  system  of  equations  is  asymptotically  stable  at  that  rest 
point. 


2.   A  method  for  computation  of  R  * 

mn 

We  assume  that  the  best  approximation  R  *  =  P*/Q*   exists,  and  is 

mn 

in  reduced  form.   The  basis  for  several  algorithms,  such  as  that  of  Remes 
(see,  for  example  [3])  and  Maehly's  second  method  [9]  is  the  following 
characterization  theorem,  which  can  be  found  in  may  texts,  e.g.  Cheney  [2]. 


Theorem  1:   In  order  that  the  irreducible  rational  function  P/Q  be  a 

best  approximation  to   f  of  the  form  R   ,   it  is  necessary  and  sufficient 
rr  mn  J 

that  the  error  have  a  least  2  +  max(m  +  deg  0,  n  +  deg  P)   alternations. 

As  with  Maehly's  second  algorithm,  we  assume  that  the  error  curve  is 

"standard"  in  that  it  has  exactly  K  =  m  +  n  +  1  points  of  zero  error 

and  K  +  1  alternations.   If  the  points  of  zero  error  are  known,  the 

best  approximation  can  be  computed  as  the  rational  function  of  the  form 

R    which  interpolates  to  the  value  of   f  at  those  points, 
mn 

Let   z  *  <  z  *  <  ...  <  z  *  be  the  points  of  zero  error  for 
1      L  K 

f  (x)  -  R  *  (x) 
E*(x)  =  -      mn 


w(x) 
Let  Z  =  (z, z„)   be  an  approximation  to   Z*  =  (z  *,..., z  *)  .   Let 

i        K  IK 

R   (x)   be  the  rational  function  which  interpolates  to   f  at   z  , ...,z   . 

mn  l      k 

Let   z  =  a,  z ■  -  =  b,   and  define  N,  =    sup     |e(x)|  ,  where  E(x)  = 
U        Ktx  k        ,     . 

zk-l-x-zk 

Rmn  ^   .   If  Z   is  close  to  Z*  ,   R    exists  and  the  N.   are 

__ mn  k 

w(x) 
finite.   Let  x,   be  the  point   in   [z,_..,z,]   such  that   |E(x  )|  =  N   , 

k  =  1,...,K  +  1  .   For  Z   close  to   Z*   these  points  must  be  unique. 


Any  of  the  variables  with  superscript  asterisks  will  denote  that  variable 
for  the  best  approximation.   Note  that  N ,  , *  -  N  *  =  0  ,  k  =  1,2,...,K  . 

K."T*X        K. 

Under  our  assumptions,  if  N   -^  ~  N,  =  0  >  k  =  1,2,...,K  ,  then  (at  least  for 

Z   in  some  neighborhood  of   Z*)   Z  =  Z*,  although  in  general  this  is  not  true. 

In  a  talk  given  at  the  joint  SIAM-MAA  meeting  for  the  Northern  California 

sections  in  February  1972,  Dr.  Milton  W.  Green  of  the  Stanford  Research 

Institute  discussed  an  algorithm  based  on  the  above  ideas  [7] .   Given  an 

initial  guess  Z  at  Z*  ,  one  computes  R    ,   and  then  N.  ,...,N..in  . 

mn  1      K+l 

The  value  of  each   z   is  then  corrected  by  changing  it  so  that  its  new 
value  is  nearer  to  x   ,   or  x   ,   as  N,    -  N   is  positive  or  negative, 
respectively.   That  is,   z,   is  changed  toward  the  point  of  largest  (in 
magnitude)  error  in  the  interval  [z,  , ,  z    ].   Dr.  Green  reported  that  he 
had  had  good  success  with  the  algorithm. 

In  an  attempt  to  systematize  the  method  and  to  make  it  amenable  to 
analysis  for  its  convergence  properties,  we  considered  the  basic  idea  in 
the  following  form.   We  formulated  the  method  as  a  continuous  (in  the  corrections 
to  the   z  )   rather  than  a  discrete  problem.   Consider  the  system  of  ordinary 
differential  equations, 

*k-"k*i-"k  >   k  =  1 K  • 

where   z,   and  N   are  as  defined  previously,  and  Z  =  Z     at   t  =  0   is 
an  approximation  to  Z*  ,  used  as  the  initial  condition.   This  system,  when 
solved  by  Euler's  method,  yields  an  algorithm  similar  to  that  proposed  by 
Dr.  Green. 


The  algorithm  we  study  is  based  on  a  slightly  different  system  of 
equations.   Although  the  convergence  properties  are  similar,  we  wished  to 
remove  the  effect  of  linear  transformations,  and  to  incorporate  some 
indication  when  the   z  *  bunch  together,  as  happens  when  f  has  a  large 
slope  at  some  point.   Consequently,  we  considered  two  somewhat  modified 
systems  of  equations,  neither  of  which  seems  to  yield  results  markedly 
superior  to  the  other.   The  first,  and  the  one  we  analyze,  is  the  system 

(1)  \  -      B       (\+i  -  V  •   k  =  1 K 

N 

where  x,   is  as  defined  previously,  and  N  =   max   N,  ,  again  with 

l<k<K+l  * 

Z  =  Z    at   t  =  0  .   The  other  system  was 

Nk+1   "  Nk    ,  , 

(2)  zk  =  2 \\  ~  ykl     i      k  =   1,...,K 

N 


where     y. 
k 


"  "fcfl  -  \ s  °  • 


fcfl  -f  "fcfl  -  \ >  °  • 


and   again     Z=Z  at      t=0. 

The  factors  xVi-i  ~  xi  »   ^nd   |  z,  -  y,  |   are  both  an  attempt  to  "slow" 
z,   when  the  z*     bunch  together.   In  some  cases  (2)  is  superior  to  (1), 
and  in  others  (1)  is  superior  to  (2) .   We  choose  to  analyze  (1)  because  it 
seems  to  be  more  consistent  in  the  optimum  "time  step"  when  solved  by  Euler's 
method.   The  analysis  of  (2)  is  nearly  the  same. 

The  point  Z*   is  clearly  a  rest  point  of  the  system  (1) ,  and  in  the 
next  section  we  give  an  analysis  showing  that  under  appropriate  assumptions, 
(1)  is  asymptotically  stable  at  Z  =  Z*  . 


There  are  many  algorithms  now  possible,  depending  on  the  numerical  method 
used  to  solve  (1).   Any  method  could  be  used,  subject  only  to  the  appropriate 
choice  of  "time  step".   However,  we  should  bear  in  mind  that  the  goal  is  not 
necessarily  to  solve  (1)  accurately,  but  rather  to  approach  Z*   closely. 
Thus,  the  use  of  Euler's  method  for  the  solution.   The  choice  of  "time  step" 
appears  to  be  a  matter  of  experience.   It  is  desired  to  use  a  near  optimal 
"time  step",  one  which  yields   R  *   to  the  desired  accuracy  in  a  minimum 
number  of  time  steps,  or  iterations.   It  is  seen  that  the  "time  step"  is  a 
parameter  similar  to  the  parameter  in  the  solution  of  elliptic  boundary 
value  problems  by  the  alternating-direction  implicit  method  [10] .   It  is 
doubtful  that  it  can  be  taken  arbitrarily  large  in  our  case,  however. 
3.   Convergence 

The  convergence  of  the  algorithms  possible  in  the  setting  of  Section  2, 
when  an  appropriate  "time  step"  is  used,  is  determined  by  whether  or  not 
the  system  (1)  is  asymptotically  stable  at  Z  =  Z*  .   We  will  make  use  of 
the  following  theorem,  which  is  paraphrased  slightly  from  the  way  it  appears 
in  most  references,  e.g.  [1], 

Theorem  _2:   The  system  of  equations 

(3)  Z  =  A  •  (Z  -  Z*)  +  H(Z  -  Z*)  , 

where  A  is  a  constant  K  x  K  matrix  and  H(Z  -  Z*)   is  a  vector  function 
which  is  small  compared  to   Z  -  Z*  ,  is  uniformly  asymptotically  stable  at  the 
point   Z  =  Z*   if  the  eigenvalues  of   A  all  have  negative  real  parts. 


Thus,  in  order  to  analyze  the  system  (1)  we  must  put  it  in  the  form  (3). 
Denote  sign  (E(x))   on   (zu_i  >  O   bY  °"k  •   By  our  assumptions  about  E(x) 


we  can  write 


K 
E(x)  =  G(x)    n   (x  -  z.)  , 
i=l       1 


where  G(x)   is  continuous  and  single  signed  if  Z   is  close  to  Z*  .   Then 

K 
we  have  N  =  a     G(x  )   f]   (xk  -  z.)  »   k=l,...,K+l.   Recall  also  that 

i=l 

a  =  zQ  <  x±   <  z±  <    .  .  .  <  xK  <  zK  <  xx+1  =  b  . 

We  now  make  an  assumption  about  the  dependence  of  G(x)   on  the  z.  . 
As  was  done  by  Maehly  and  Witzgall  [9],  we  assume  that  near  the  point  Z*  , 
the  function  G(x)   does  not  depend  very  much  on  the   z.  .   (Note:   In  the 
case  of  approximating  x   by  a  polynomial  of  degree  <  K  -  1  ,   with 
w(x)  =  1  ,  we  have  G(x)  =  1.)   Then  we  have,  for  a  given  value  of  x  , 

9E(x)     n/  :   k  ,  N      E(x) 

zt~  a  "  G(x)  ,"   (x  -  zi}  =  -  irrr.    • 


Then,  since  x   is  the  point  of  extreme  error,  even  though  x,   may  change 
significantly  with   z.  ,   it  is  seen  that  the  extreme  value  does  not.   Thus 

9Nk  °k   E(xk}        N 

we  have  - —  ~ ~  ,   since  near  Z*  ,   N.  «  N.   Then  we 

9z,  x.-z.    x.  -  z.                     k 

j  k    j     k    j 

have 


+  Hk  ,  k  =  1,...,K  , 


Hl 

where  H  =  ^  !  J      is  small  compared  to   Z  -  Z*  .   Further  simplification  gives 
HK 


the  expression 


K  * 

z  .  -  z  * 


(x.  ^_  *  -  x  *)  /      -. -1  ■} — — -  +  H,  ,  k  =  1 K  . 

k+1     k    JL^   (x.  _*  -  z.*)  (x  *  -  z.*)  k  '       '    ' 
•  i    k+1     j     k     j 
J=l 


Thus,  we  have  the  system  (1)   in  the  form  (3)  with 

2 

Cx   *  -  x  *) 

^xi,+i     xi,  i 

(4)         A  = 


(x    *   -   z   *)    (x   *    -   z   *) 

K  k+1        j  ;  ^  k        j  ; 


k,j  =  1.....K  . 


Now  we  must  investigate  the  location  of  the  eigenvalues  of  the  matrix 
A  .   To  simplify  the  notation,  we  drop  the  asterisks  from  the  variables.   We 
will  show  that  the  matrix  -A   is  a  matrix  of  class  |(   (see  Fiedler  &  Ptak 
[6].   These  matrices  are  also  known  as  M-matrices.) 

Let  B  be  a  square  matrix  with  nonpositive  elements,  except  possibly 
on  the  diagonal.   Matrices  of  class   ((  are  a  subset  of  such  matrices, 
characterized  by  many  equivalent  properties  [6].   The  two  equivalent  pro- 
perties we  shall  use  are  given  in 

Theorem  3:  The  following  properties  of  B  are  equivalent.  (i)  The  real 
part  of  every  eigenvalue  of  B  is  positive,  (ii)  Every  leading  principal 
minor  of   B   is  positive. 

Lemma  4:   The  matrix  A  defined  by  (4)  has  negative  elements  on  the  diagonal 
and  positive  elements  elsewhere.   Further,  the  principal  minors  of   -A  , 


(xk+l  "  xk} 

Det  |-  (     _  z  )  (   _  Z)J  are  positive  for   f-  1.2.....K  , 

whenever  X;L  <  Z;L  <  x2  <  z2  <  .  .  .  <  xR  <  zR  <  xR+1  . 

(xk+l  "  Xk} 

Proof:   We  consider  the  element  a.  .  =  -, : — ; -  .   The  numerator 

kj    (xk+1  -  z.)  (xk  -  z.) 

is  positive,  and  if   j  <  k  ,   then   z .  <  x .  .  n  <  x.  <  x.  . .  ,   and  hence 

J     J+l  ~  k    k+1 

a.  .  >  0  .   If   j>k,we  have  x,  <  x.  ,  n  <  x .  <  z .  ,   and  again  a,  .  >  0  . 
kj  k    k+1  -   j    j  kj 

If   i  =  k  ,  x.  =x.  <  z.  <  x,,.  =  x.  ,,  ,   hence  a,  .  <  0  .   We  now  consider 
k    j    j    j+l    k+1  kj 

i     „th  -  . 

the  I  principal  minor  of   -A  , 

(xk+l  "  Xk)2 
Det 


(xk+l  "  Zi}  (xk  ~  z1} 

1    J    k    J  /   k.j  -  1,...,* 

2 
Removal  of  the  factor   (x,,,  -  x,  )    from  each  row  does  not  change  the  sign 

of  the  determinant,  so  we  consider  the 


I  k.j  =  l,...,jg 


The  value  of  the  latter  determinant  being  nonzero  is  equivalent  to  the 
existence  of  a  unique  solution  to  the  following  interpolation  problem:   With 
functions  gk(y)  =  -  7^ x  , — — r  ,  k  =  l,...,i   ,   and  given  points 

(y.»w.)  ,   j  =  1.....J0  ,   find  constants  o>       k  =  l.....i   such  that 
J    J  k 


'\+l   'j'  v"k   'j 


^k 


W  =  §"  <x^  -  yj  (x.  -  y J  =  wj  * 


Here  we  assume  that  none  of  the  y.   coincide  with  any  of  the  x,  ,  which 

3  k 

we  have  guaranteed  in  the  case  of  interest.   The  above  discussion  says  that 

D    is  nonzero  if  the  set  of  functions   g,  (y) ,  k  =  1 , . . . , I    is  unisolvent  on 

the  permitted  set  of  points.   See  Davis  [5]  for  further  discussion. 

The  interpolation  problem  is  known  to  be  uniquely  solvable,  if  and  only 

if  any  linear  combination  of  the   g^(y)  >  k  =  1,...,£  which  is  zero  at  I 

distinct  points  is  identically  zero.   Consider  G  (y)  =  /J  ai_8i,(y)  •   Now 

k=l 

i+i  P„  ,  (y) 

c0(y)  -  - 


1 Z  ak   n        <x,  -  y)  ■ 


. n,  (xk  -  y)  ^k 

k=1  U+i 


where  P„_,   and  Q     are  polynomials  of  degree   £  -  1   and  I  +   1  ,   res- 
pectively.   We  see  that  P   1  ,   and  hence  G   can  have  at  most  I   -   1     dis- 
tinct zeros,  unless  G   is  identically  zero.   Hence  D   *  0  . 

Now  let  y  ,...,y   be  variable,  but  such  that  each  y.   satisfies  the 
1      *  J 

condition  stated  for   z.   in  the  lemma,  i.e.,   x.  <  y.  <  x.in.   Let  y. y, 

J  j    j    j+1         1      I 

replace   z  ,...,z   ,   respectively  in  D   ,   and  note  that  x  -  y,   appears 

X  X/  X-  K.       K. 

only  on  the  diagonal.   Hence  by  choosing  y,   sufficiently  close  to  x,  ,  the 
determinant  is  diagonally  dominant,  and  is  thus  positive.   By  the  continuity 
of  D   as  a  function  of   z  ,...,z   ,   and  the  fact  D   is  never  zero,  we 
conclude  that   D   >  0  .   This  completes  the  proof  of  the  lemma. 

Application  of  Theorems  3  and  2  to  the  system  (1)  yields  the  fact  that 
(1)  is  asymptotically  stable  at  Z  =  Z*  . 


4.   Numerical  Implementation. 

A  version  of  the  algorithm,  which  we  will  call  "Algorithm  G",  was 

implemented  in  double  precision  Fortran  on  the  IBM  360/67  at  the  Naval 

Postgraduate  School.   We  used  Euler's  Method  to  solve  (1)  numerically,  using 

the  initial  guess  Z     to  be  the  zeros  of  the  Chebyshev  polynomial  of 

degree  K  ,  translated  to  the  interval   [a,b] .   The  rational  function  R 

mn 

was  found  by  solving  the  linear  system  obtained  by  requiring  R    to 

interpolate  to   f  at  the   z.  .   The  IMSL   routine  LEQTIF  was  used  to 

solve  the  system. 

The  extreme  values  N,   were  found  by  the  method  suggested  by  Maehly 

and  Witzgall  [9].   A  search  is  made  for  a  "turning  point"  using  the 

previous  value  of  x,   as  an  initial  estimate,  using  steps  of 

h  =  .  015(z,  -  z.  -).   When  three  points  have  been  computed  so  that 
k         k    k-1 

|E(x)|   is  largest  at  the  middle  point,  the  value  of  x^  is  approximated 

by  passing  a  parabola  through  the  three  points,  and  finding  its  extreme 

point.   The  alternative  is  to  convert  the  problem  to  a  discrete  one  by 

evaluating  the  error  at  a  fixed  number  of  points,  as  did  Lee  and  Roberts  [8] 

We  feel  our  method  is  probably  faster,  especially  in  the  latter  stages, 

although  it  assumes  the  error  curve  has  but  one  local  extreme  in  each 

interval   (zi_i  >  ziJ  •  ^e  believe  the  method  to  be  more  accurate  for  smooth 

problems,  as  well  as  preserving  the  continuity  of  the  original  problem. 


1  IMSL  -  International  Mathematical  and  Statistical  Libraries,  Inc., 
6200  Hillcroft,  Suite  510,  Houston,  Texas   77036 
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Algorithm  G  was  tested  by  running  a  variety  of  problems.   The  same 

2 
problems  were  also  run  using  the  Remes  algorithm  and  Maehly's  second 

method;  the  latter  program  was  written  by  the  author.   The  problems  ranged 

from  "easy",  such  as  exp(x),   T(x) ,  and  log  x,  to  "hard"  problems  such 

I —  /—      1         1 

as  vx  on  [0,1]  and  vx  on  ["Tjl]  and  [y-r,l]  ,  the  latter  two  with  weight 

function  w(x)  =  Vx-  for  relative  fits.   In  addition,  the  function  r(x)  = 
(arctan  8x)   v(8x-l)   +  1  /8x  ,   with  w(x)  =  1  ,  on  [-1,1]  was  attempted. 
This  latter  function  is  due  to  Rutishauser  (see  Cody  &  Stoer  [4],  p.  179) 
and  is  difficult  because  the  initial  guesses  usually  used  (data  from 
appropriate  Chebyshev  polynomial)  with  the  Remes  algorithm  lead  to 
approximants  with  a  pole  in  [-1,1] . 

The  approximations  of  the  form  R.  ,  R  _ ,  R?9,  and  R,„  were  usually 

attempted,  although  in  some  specific  instances  others  were  computed.   In 
the  case  of  the  "easy"  functions  all  three  algorithms  worked  well,  with 
the  Remes  algorithm  converging  very  rapidly,  of  course.   The  Remes  algorithm 
and  Maehley's  method  failed  on  one  or  more  of  the  "hard"  problems,  while 
Algorithm  G,  for  appropriate  "time  step",  did  not  fail.   We  note  that  one 
iteration  of  the  Remes  algorithm  requires  the  solution  of  a  non-linear 
system,  Maehly's  method  requires  the  solution  of  two  linear  systems,  while 
Algorithm  G  requires  the  solution  of  one  linear  system. 

The  difficulties  with  Rutishauser 's  function  r(x)   were  investigated 
more  closely.   Aside  from  the  problem  of  poles,  the  function  has  a  large 
slope  near  x  =  0  ,  which  apparently  causes  difficulties.   The  programs 


2  The  IMSL  routine  IRATCU  was  used.   This  is  a  Fortran  version  of  procedure 
Che.by6h.ev   due  to  Cody,  Fraser,  and  Hart  [3]. 
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for  the  Remes  algorithm  and  Maehly's  method  were  modified  to  accept  input 
initial  guesses.   For  initial  guesses  at  the  extremes  which  were  accurate 
to  seven  significant  digits,  the  Remes  algorithm  failed.   For  initial  guesses 
at  the  interpolation  points  which  were  accurate  to  seven  significant  digits, 
Maehly's  method  failed.   No  particular  difficulties  were  experienced  by 
Algorithm  G. 

With  regard  to  the  possibility  of  poles  in  the  initial  approximant, 
we  have  discovered  that  while  theoretically  the  method  fails,  the  numerical 
algorithm  may  be  able  to  recover.   For  example,  when  approximating  sin  x 
on  [0,  4.1]  with  R   =  P_/Q   ,  using  the  Chebyshev  points  as  initial 

guesses  gives  an  initial  approximant  with  a  pole  near  x  =  1.7.   However, 
because  of  the  approximation  to  the  extremes,  the  routine  recovers,  forces 
the  pole  out  of  the  approximation,  and  then  converges  to  the  correct  result. 
This  may  not  be  a  general  rule,  however,  but  is  an  interesting  example  of 
how  robust  the  algorithm  can  be. 

Having  satisfied  ourselves  that  the  algorithm  works  quite  well,  a  study 
was  made  of  how  one  should  choose  the  "time  step".   We  found  that  "time  step" 
At  =  .20  was  usually  (not  always)  small  enough  for  convergence.   The  optimum 
value  for  At   is  dependent  on  the  problem,  and  sometimes  is  significantly 
larger  than   .20.   As  might  be  expected,  the  optimum  value  is  larger  than 
required  for  accurate  solution  of  (1) ,  and  it  is  better  to  underestimate 
At   than  to  overestimate  it.   Table  1  gives  the  number  of  time  steps  (or 
iterations)  required  for  convergence  of  various  approximations  to  exp(x) 
on  [0,1],   Vx~  on  [0,1],  and  r(x)   on  [-1,1],  versus  At  .   The  iterations 
were  stopped  when  successive  approximations  were  obtained  whose  respective 
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coefficients  had  relative  differences  of  less  than  10 

Because  the  solution  of  (1)  is  close  to  the  solution  of  Z  =  A(Z  -  Z*) 
for  large  times,  the  convergence  rate  is  seen  to  be  linear.   Further,  we 
can  see  in  Table  1  that  convergence  is  slow  even  for  smooth  functions 
such  as  exp(x),  where  10-20  iterations  are  required.   In  the  case  of  v^~  > 
where  most  of  the  z   tend  to  bunch  near  zero,  significantly  more  iterations 
are  required,  the  number  increasing  with  K  . 

Table  1:   Search  for  optimum  At 
exp(x)        \/x~ r(x) 


\t 

VQ.1 

VQ2 

VQ1 

P2/Q2 

VQ.1 

P2/Q2 

VQ2 

.10 

36 

64 

.125 

51 

>100 

29 

64 

61 

.15 

23 

47 

41 

97 

22 

49 

53 

.175 

44 

82 

17 

39 

42 

.20 

16 

35 

>100 

>100 

22 

failed 

failei 

.225 

14 

32 

.25 

12 

28 

.275 

19 

.30 

36 

23 

.325 

.35 

19 

.375 

18 

.40 

16 

.425 

16 
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5.  Conclusion. 

Algorithm  G  appears  to  be  very  robust.   When  coupled  with  the  appropriate 
"time  step",  it  is  likely  infallible,  provided  the  initial  guess  did  not 
yield  an  approximant  with  a  pole  in  the  interval.   Even  then,  the  version 
of  the  algorithm  we  implemented  has  recovered  in  specific  instances. 

While  this  algorithm  cannot  compete  with  the  Remes  algorithm  on  the 
basis  of  speed,  our  tests  show  its  speed  to  be  comparable  to  Maehly's 
algorithm,  which  shows  up  quite  well,  on  that  basis,  in  the  Lee  and  Roberts 
study.   One  suspects  the  Lee  and  Roberts  timing  is  biased  since  cases  where 
the  algorithm  was  successful  were  likely  to  be  the  relatively  easier  (and 
faster)  cases. 

The  simplicity  of  the  algorithm  coupled  with  its  high  success  rate, 
make  it  worthy  of  consideration  for  computing  approximations  which  result 
in  failure  of  the  Remes  algorithm.   The  principal  disadvantage  seems  to 
be  the  assumption  that  there  are  no  more  than  m  +  n  +  1   zeros  of  the 
error  curve.   Degenerate  cases  can  be  handled  by  increasing  the  degree  of 
the  numerator  and/or  denominator  so  that  the  error  curve  becomes  standard. 
If  the  cause  for  more  than  m  +  n  +  1  zeros  is  not  degeneracy,  the  algorithm 
may  fail. 
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